<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
    <title>bvodeS</title>
  </head>
  <body bgcolor="#FFFFFF">
    <center>Scilab Function</center>
    <div align="right">Last update : 08/02/2006</div>
    <p>
      <b>bvodeS</b> - simplified call of bvode</p>
    <h3>
      <font color="blue">Calling Sequence</font>
    </h3>
    <dl>
      <dd>
        <tt>z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,
         [ystart,dfsub,dgsub,fixpnt,ndimf,ndimi,ltol,tol,ntol,nonlin,
          collpnt,subint,iprint,ireg,ifail])
    </tt>
      </dd>
    </dl>
    <h3>
      <font color="blue">Parameters</font>
    </h3>
    <ul>
      <li>
        <tt>
          <b>x</b>
        </tt>:  array of points at which approximations of the solution will be computed.
    The points x(i) must be given in increasing order.</li>
      <li>
        <tt>
          <b>m</b>
        </tt>: array of orders for the differential equations given in <tt>
          <b>fsub</b>
        </tt>.</li>
      <li>
        <tt>
          <b>n</b>
        </tt>:  number of differential equations, length of <tt>
          <b>m</b>
        </tt>.</li>
      <li>
        <tt>
          <b>a</b>
        </tt>: left end of solution interval, where a&lt;=x(1).</li>
      <li>
        <tt>
          <b>b</b>
        </tt>:  right end of solution interval, where x($)&lt;=b.</li>
      <li>
        <tt>
          <b>fsub</b>
        </tt>:  name of a Scilab function for evaluating the rhs of the differential 
       equation system to be solved; fsub also may be a list for parameter transfer.</li>
      <li>
        <tt>
          <b>gsub</b>
        </tt>:   name of a Scilab function for evaluating the boundary conditions for the
       given differential equation system; gsub may also be a list as for fsub.</li>
      <li>
        <tt>
          <b>zeta</b>
        </tt>: array of points at which the boundary or side conditions are given in 
       increasing order. Each point must occur as often as boundary or side 
       conditions are given there. Each side condition point other than a or b must
       be included in the array fixpnt.
       The length of zeta should be sum(m).</li>
      <li>
        <tt>
          <b>z</b>
        </tt>:  array of all derivatives of the solution functions up to the order
    m(i)-1 for the i-th function. (See the examples below.) The length of z should
    be sum(m).</li>
      <li>
        <tt>
          <b>ystart,dfsub,dgsub,fixpnt,ndimf,ndimi,ltol,tol,ntol,nonlin,
          collpnt,subint,iprint,ireg,ifail</b>
        </tt>: These optional arguments may be called by name in any order in the
form argument=name. The meaning of ystart to ireg is given in the <a href="bvode.htm">
          <tt>
            <b>bvode</b>
          </tt>
        </a> help page.
The Scilab functions <tt>
          <b>dfsub</b>
        </tt> and <tt>
          <b>dgsub</b>
        </tt> for evaluating the Jacobians may also be lists
for parameter transfer. The function <tt>
          <b>ystart</b>
        </tt> is called <tt>
          <b>guess</b>
        </tt> in bvode.</li>
      <li>
        <tt>
          <b>ifail</b>
        </tt>:  if ifail=1, all parameters needed for the call of bvode are displayed.</li>
    </ul>
    <h3>
      <font color="blue">Description</font>
    </h3>
    <p>
This interface program simplifies the call of <a href="bvode.htm">
        <tt>
          <b>bvode</b>
        </tt>
      </a>, a program for the numerical
solution of multi-point boundary value problems for mixed order systems of ordinary
differential equations. The Scilab program bvode is adapted from the fortran program
colnew. See the paper by U. Asher, J. Christiansen, and R.D. Russel:
Collocation software for boundary-value ODE's, ACM Trans. Math. Soft. 7:209-222, 1981.
The following examples should demonstrate not only how such systems can be solved 
with the help of bvodeS, but also should emphazise some important problems which occur 
with boundary value problems for ordinary ode's.
</p>
    <h3>
      <font color="blue">Examples</font>
    </h3>
    <pre>
// 1. Modified example from help bvode. 

// DE:  y1''''(x)=(1-6*x*x*y1'''(x)-6*x*y1''(x))/(y2(x)^3)
//         y2'(x)=1
// z=[y1 ; y1' ; y1'' ; y1''' ; y2]
// BV: y1(1)=0 y1''(1)=0 
// BV: y1(2)=1 y1''(2)=0 y2(2)=2

function RhS=fsub(x,z)
  RhS=[(1-6*x*x*z(4)-6*x*z(3))/(z(5)^3);1]
endfunction

function g=gsub(i,z)
  g=[z(1) z(3) z(1)-1 z(3) z(5)-2]
  g=g(i)
endfunction

function [z,lhS]=ystart(x)
  z=zeros(5,1);z(5)=1;
  lhS=[0;1];
endfunction

n=2;
m=[4 1];
N=100;
a=1; b=2;
zeta=[a a b b b];
x=linspace(a,b,N);

ltol=4; // We want to change the default error for y1'''.
tol=1e-12;
tic()
z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ltol=ltol,tol=tol,ystart=ystart);
// Try tol=1e-14 etc.
toc()

function z=yex(x) // True solution
  z=zeros(5,1);
  z(1)=0.25*(10*log(2)-3)*(1-x)+0.5*(1/x+(3+x)*log(x)-x)+(x-1)
  z(2)=-0.25*(10*log(2)-3)+0.5*(-1/x^2+(3+x)/x+log(x)-1)+1
  z(3)=0.5*(2/x^3+1/x-3/x^2)
  z(4)=0.5*(-6/x^4-1/x/x+6/x^3)
  z(5)=x
endfunction

zex=[];for xx=x, zex=[zex yex(xx)]; end
scf(0); clf();
plot2d(x,abs(z-zex)',style=[1 2 3 5 6])
xtitle('Absolute error','x',' ')
legend(['z1(x)';'z2(x)';'z3(x)';'z4(x)';'z5(x)'])

// example #2. An eigenvalue problem

// y''(x)=-la*y(x)
// BV: y(0)=y'(0); y(1)=0
// Eigenfunctions and eigenvalues are y(x,n)=sin(s(n)*(1-x)), la(n)=s(n)^2,
// where s(n) are the zeros of f(s,n)=s+atan(s)-(n+1)*pi, n=0,1,2,...
// To get a third boundary condition, we choose y(0)=1
// (With y(x) also c*y(x) is a solution for each constant c.)
// We solve the following ode system:
// y''=-la*y
// la'=0
// BV: y(0)=y'(0), y(0)=1; y(1)=0
// z=[y(x) ; y'(x) ; la]

function rhs=fsub(x,z)
  rhs=[-z(3)*z(1);0]
endfunction

function g=gsub(i,z)
  g=[z(1)-z(2) z(1)-1 z(1)]
  g=g(i)
endfunction

// The following start function is good for the first 8 eigenfunctions.
function [z,lhs]=ystart(x,z,la0)
  z=[1;0;la0]
  lhs=[0;0]
endfunction

a=0;b=1;
m=[2;1];
n=2;
zeta=[a a b];
N=101;
x=linspace(a,b,N)';

// We have s(n)-(n+1/2)*pi -&gt; 0 for n to infinity.
la0=input('n-th eigenvalue: n= ?');la0=(%pi/2+la0*%pi)^2;

z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,la0));

clf()
plot2d(x,[z(1,:)' z(2,:)'],style=[5 1],axesflag=5) 
xtitle(['Startvalue =  '+string(la0);'Eigenvalue = '+string(z(3,1))],'x',' ')
legend(['y(x)';'y''(x)'])


// example #3. A boundary value problem with more than one solution.

// DE: y''(x)=-exp(y(x))
// BV: y(0)=0; y(1)=0
// This boundary value problem has more than one solution.
// It is demonstrated how to find two of them with the help of
// some preinformation of the solutions y(x) to build the function ystart.
// z=[y(x);y'(x)]

a=0;b=1;m=2;n=1;
zeta=[a b];
N=101;
tol=1e-8*[1 1];
x=linspace(a,b,N);

function rhs=fsub(x,z),rhs=-exp(z(1));endfunction

function g=gsub(i,z)
  g=[z(1) z(1)]
  g=g(i)
endfunction

function [z,lhs]=ystart(x,z,M) 
  //z=[4*x*(1-x)*M ; 4*(1-2*x)*M]
  z=[M;0]
  //lhs=[-exp(4*x*(1-x)*M)]
  lhs=0
endfunction

for M=[1 4]
   if M==1
      z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,M),tol=tol);
   else
      z1=bvodeS(x,m,n,a,b,fsub,gsub,zeta,ystart=list(ystart,M),tol=tol);
   end
end

// Integrating the ode yield e.g. the two solutions yex and yex1. 

function y=f(c),y=c.*(1-tanh(sqrt(c)/4).^2)-2;endfunction 
c=fsolve(2,f);

function y=yex(x,c)
  y=log(c/2*(1-tanh(sqrt(c)*(1/4-x/2)).^2))
endfunction 

function y=f1(c1), y=2*c1^2+tanh(1/4/c1)^2-1;endfunction
c1=fsolve(0.1,f1);

function y=yex1(x,c1)
  y=log((1-tanh((2*x-1)/4/c1).^2)/2/c1/c1)
endfunction 

disp(norm(z(1,:)-yex(x)),'norm(yex(x)-z(1,:))= ')
disp(norm(z1(1,:)-yex1(x)),'norm(yex1(x)-z1(1,:))= ')

clf();
subplot(2,1,1)
plot2d(x,z(1,:),style=[5])
xtitle('Two different solutions','x',' ') 
subplot(2,1,2)
plot2d(x,z1(1,:),style=[5])
xtitle(' ','x',' ')


// example #4. A multi-point boundary value problem.

// DE y'''(x)=1
// z=[y(x);y'(x);y''(x)]
// BV: y(-1)=2 y(1)=2
// Side condition: y(0)=1 

a=-1;b=1;c=0;
// The side condition point c must be included in the array fixpnt.
n=1;
m=[3];

function rhs=fsub(x,z)
  rhs=1
endfunction

function g=gsub(i,z)
  g=[z(1)-2 z(1)-1 z(1)-2]
  g=g(i)
endfunction

N=10;
zeta=[a c b];
x=linspace(a,b,N);

z=bvodeS(x,m,n,a,b,fsub,gsub,zeta,fixpnt=c);
          
function y=yex(x)
y=x.^3/6+x.^2-x./6+1
endfunction

disp(norm(yex(x)-z(1,:)),'norm(yex(x)-z(1,:))= ')
 
  </pre>
    <h3>
      <font color="blue">See Also</font>
    </h3>
    <p>
      <a href="bvode.htm">
        <tt>
          <b>bvode</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="ode.htm">
        <tt>
          <b>ode</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="dassl.htm">
        <tt>
          <b>dassl</b>
        </tt>
      </a>,&nbsp;&nbsp;</p>
    <h3>
      <font color="blue">Author</font>
    </h3>
    <p>Rainer von Seggern</p>
  </body>
</html>
